package org.farrell.EOSCalc;

import static java.lang.Math.PI;
import static java.lang.Math.abs;
import static java.lang.Math.exp;
import static java.lang.Math.log10;
import static java.lang.Math.pow;
import static java.lang.Math.signum;
import static java.lang.Math.sqrt;
import java.util.Arrays;
import org.apache.commons.math.complex.Complex;

import android.os.Bundle;

public abstract class CubicEOS {
	// Implements a generic Cubic Equation of State for a single compound
	// internally the class uses SI units: Temperature in Kelvin, Pressure in Pa
	// Energy in J, etc.
	
	// Constants and properties of the compound
	public static final double R = 8.3143; //J/mol*K
	public final double Tc, Pc;
	public final double w;
	
	// Errors and iteration levels for saturation calculations
	public static final int MAX_ITER = 100;
	public static final double MAX_ERR = 1e-12;
	
	// Acceptable deviation of the gibbs' departure function in the 2-root region
	public static final double VLE_REGION = 0.01;
	
	// The keys to the bundle
	public static final String KEY_NAME = "name";
	public static final String KEY_TEMP = "T";
	public static final String KEY_PRES = "P";
	public static final String KEY_ROOTS = "roots";
	public static final String KEY_Z = "Z";
	public static final String KEY_VOL = "V";
	public static final String KEY_FUG = "f";
	public static final String KEY_GIBBS = "G";
	public static final String KEY_ENTHALPY = "H";
	public static final String KEY_ENTROPY = "S";
	
	public CubicEOS(double criticalTemperature, double criticalPressure, double acentricFactor) {
		// Initialise the compound, Temp is in Kelvin, Pressure in Pascals
		Tc = criticalTemperature;
		Pc = criticalPressure;
		w = acentricFactor;
	}
	// The Abstract methods that define a particular cubic eos
	public abstract double getA(double Temperature, double Pressure);
	
	public abstract double getB(double Temperature, double Pressure);
	
	public abstract Double[] getZ(double a, double b);
	
    public abstract double getHdepart(double z, double a, double b, double temp);
    
    public abstract double getSdepart(double z, double a, double b, double temp);
    
    public abstract double getGdepart(double z, double a, double b, double temp);
    
    
	// The common methods to all cubic eos
	public Bundle setState(double Temperature, double Pressure){
		// Given a Temperature and Pressure calculates all the component properties
		// returns a bundle with the appropriate arrays loaded in.
		double A = getA(Temperature, Pressure);
		double B = getB(Temperature, Pressure);
		Double[] Zinterim = getZ(A, B);
		Double[] Z = findRoots(Zinterim, A, B, Temperature);
		Bundle results = new Bundle();
		int numRoots = Z.length;
		double[] z = new double[numRoots];
		double[] v = new double[numRoots];
		double[] f = new double[numRoots];
		double[] g = new double[numRoots];
		double[] h = new double[numRoots];
		double[] s = new double[numRoots];
		for (int i=0; i<numRoots; i++) {
			z[i] = Z[i].doubleValue();
			v[i] = getV(z[i],Temperature,Pressure);
			g[i] = getGdepart(z[i], A, B, Temperature);
			f[i] = getFugacity(g[i], Pressure);
			g[i] = R*Temperature*g[i];
			h[i] = R*Temperature*getHdepart(z[i], A, B, Temperature);
			s[i] = R*getSdepart(z[i], A, B, Temperature);
		}
		results.putDouble(KEY_TEMP, Temperature);
		results.putDouble(KEY_PRES, Pressure);
		results.putInt(KEY_ROOTS, numRoots);
		results.putDoubleArray(KEY_Z, z);
		results.putDoubleArray(KEY_VOL, v);
		results.putDoubleArray(KEY_FUG, f);
		results.putDoubleArray(KEY_GIBBS, g);
		results.putDoubleArray(KEY_ENTHALPY, h);
		results.putDoubleArray(KEY_ENTROPY, s);
		return results;
	}
                
    public Double[] cubicRoots(double a0, double a1, double a2){
        double Q = (a2*a2-3*a1)/9.0;
        double R = (27*a0 + a2*(2*a2*a2-9*a1))/54.0;
        double t = a2/3.0;
        Double[] z;
        
        if ((R*R)<(Q*Q*Q)){
        	// There are three real roots, and thus at most 2 stable roots
        	
        	// setting up complex versions of R and Q
        	Complex cR = new Complex(R, 0.0);
        	Complex cQ = new Complex(Q, 0.0);
        	
        	// some intermediary complex numbers
        	Complex s = new Complex(2*PI/3, 0.0);
        	Complex cQ3 = cQ.multiply(cQ.multiply(cQ));
        	Complex m = cR.divide(cQ3.sqrt());
        	
        	// using trigonometry to find the roots
        	Complex theta = m.acos().multiply(1/3.0);
        	Complex x1 = cQ.sqrt().multiply(-2.0).multiply(theta.cos());
        	Complex x2 = cQ.sqrt().multiply(-2.0).multiply(theta.add(s).cos());
        	Complex x3 = cQ.sqrt().multiply(-2.0).multiply(theta.subtract(s).cos());
        	
        	// getting the real results
        	double z1 = x1.getReal() - t;
        	double z2 = x2.getReal() - t;
        	double z3 = x3.getReal() - t;
        	double[] zinterim = {z1, z2, z3};

        	// sorting the roots and ditching the middle root, which is unphysical
        	Arrays.sort(zinterim);
        	z = new Double[2];
        	z[0] = zinterim[0];
        	z[1] = zinterim[2];
        } else {
        	// there is only one real root
        	z = new Double[1];
        	
        	if ((R!=0)&&(Q!=0)){
        		double A = -signum(R)*pow(abs(R)+sqrt(R*R-Q*Q*Q),1/3.0);
        		double B = Q/A;
        		z[0] = (A+B)-t;
        	} else {
        		z[0] = -t;
        	}
        }

        return z;
	}

	public double[] getLnPhi(Double[] z, double a, double b, double temperature) {
		// Calculates the log of the fugacity coefficient for each root
		// Returns an array of size 2 with the results
		double[] G = new double[z.length];
		for (int i =0; i<G.length; i++){
			G[i] = getGdepart(z[i],a,b,temperature);
		}
		return G;
	}
	
	public Double[] findRoots(Double[] zinterim, double a, double b, double temperature) {
		// Determines how many stable roots exist (1 or 2) if the cubic equation
		// returns two roots
		Double [] z = null;
		Arrays.sort(zinterim);
		if (zinterim.length>1){
			if (zinterim[1]<b){
				// At high pressures the roots are less than B, in this case 
				// only the largest root is physical
				z = new Double[1];
				z[0] = zinterim[1].doubleValue();
			} else {
				// Need to check which root is most thermodynamically stable
				double [] lnPhi = getLnPhi(zinterim, a ,b, temperature);
				if (abs(lnPhi[0]-lnPhi[1])>abs(VLE_REGION*lnPhi[0])) {
					// Within 1% difference, probably in the 2 root region
					z = new Double[1];
					if(lnPhi[0] > lnPhi[1]) {
						// Gas phase is most stable
						z[0] = zinterim[1];
					} else {
						// liquid phase is most stable
						z[0] = zinterim[0];
					}
				} else {
					return zinterim;
				}
			}
			return z;
		}
		return zinterim;
	}

	public double diff_G(double Temperature, double Pressure) throws VLEException {
		double a = getA(Temperature, Pressure);
		double b = getB(Temperature, Pressure);
		Double[] z = getZ(a,b);
		double[] g = getLnPhi(z,a,b,Temperature);
		double diff = 0.0;
		try {
			diff = g[1]-g[0];
		}
		catch (java.lang.ArrayIndexOutOfBoundsException e){
			VLEException err = new VLEException();
			err.putMessage("Failure to locate two phase region");
			throw err;
		}
		return diff;
	}

    public double guessPsat(double Temperature){
        // Helper function, uses a crude correlation from the Clapeyron equation
        // to guess at a value for Psat
        double logP = (7.0/3.0)*(w+1)*(1-(Tc/Temperature));
        double Pr = pow(10.0,logP);
        return Pr*Pc;
    }

    public double guessTsat(double Pressure){
        // Helper function, uses a crude correlation from the Clapeyron equation
        // to guess at a value for Tsat, P is in Pa
        double Pr = Pressure/Pc;
        double Tr = 1-(3.0/7)*log10(Pr)/(w+1);
        return Tc/Tr;
    }

    public double getPsat(double Temperature) throws VLEException {
        // Finds the saturation pressure using the guess function and false-position
    	// Takes temperature in K returns a pressure in Pa
        if (Temperature > Tc) {
        	VLEException error = new VLEException();
        	error.putMessage("A two phase region does not exist when the temperature is above the Critical Temperature.");
            throw error;
        }
        double guess = guessPsat(Temperature);
        double dfguess = diff_G(Temperature, guess);
        double low = 0.0, high = 0.0;
        if (dfguess < 0){
            while (dfguess<0){
                low = guess;
                guess = 1.05*guess;
                dfguess = diff_G(Temperature, guess);
            }
            high = guess;
        }
        else if (dfguess > 0){
            while (dfguess>0){
                high = guess;
                guess = 0.95*guess;
                dfguess = diff_G(Temperature, guess);
            }
            low = guess;
        }
        else {
            return guess;
        }
        
        double f_low = diff_G(Temperature,low);
        double f_high = diff_G(Temperature,high);
        double f_half = 0.0;
        double err = 10, half = 0;
        int count = 0;
        while ((err > MAX_ERR)&&(count<MAX_ITER)){
            half = f_low*(low-high)/(f_high-f_low)+low;
            f_half = diff_G(Temperature,half);
  
            if (f_half < 0){
                low = half;
                f_low = diff_G(Temperature,low);
            }
            else if (f_half > 0){
                high = half;
                f_high = diff_G(Temperature,high);
            }
            else{
                return high;
            }
            count++;
            err = abs(f_half);
        }
        return half;
    }
    
    public double getTsat(double Pressure) throws VLEException{
        // Finds the saturation temperature using the guess function and false-position
    	// Takes pressure in Pa returns temperature in K
        if (Pressure > Pc) {
        	VLEException error = new VLEException();
        	error.putMessage("A two phase region does not exist when the pressure is above the Critical Pressure.");
            throw error;
        }
        double guess = guessTsat(Pressure);
        double dfguess = diff_G(guess, Pressure);
        double low = 0.0, high = 0.0;
        if (dfguess > 0){
            while (dfguess>0){
                low = guess;
                guess = 1.05*guess;
                dfguess = diff_G(guess, Pressure);
            }
            high = guess;
        }
        else if (dfguess < 0){
            while (dfguess<0){
                high = guess;
                guess = 0.95*guess;
                dfguess = diff_G(guess, Pressure);
            }
            low = guess;
        }
        else {
            return guess;
        }
        
        double f_low = diff_G(low, Pressure);
        double f_high = diff_G(high, Pressure);
        double f_half = 0.0;
        double err = 10, half = 0;
        int count = 0;
        while ((err > MAX_ERR)&&(count<MAX_ITER)){
            half = f_low*(low-high)/(f_high-f_low)+low;
            f_half = diff_G(half, Pressure);
  
            if (f_half < 0){
                low = half;
                f_low = diff_G(low, Pressure);
            }
            else if (f_half > 0){
                high = half;
                f_high = diff_G(high, Pressure);
            }
            else{
                return high;
            }
            count++;
            err = abs(f_half);
        }
        return half;
    }

    public double getFugacity(double lnphi, double Pressure){
    	return exp(lnphi)*Pressure;
    }
    
    public double getV(double z, double Temperature, double Pressure){
    	return (z*R*Temperature)/Pressure;
    }
    
}
